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We present cosmological bounds on the thermal axion mass in an extended cosmological scenario in 
which the primordial power spectrum of scalar perturbations differs from the usual power-law shape 
predicted by the simplest inflationary models. The power spectrum is instead modeled by means 
of a ’’piecewise cubic Hermite interpolating polynomial” (PCHIP). When using Cosmic Microwave 
Background measurements combined with other cosmological data sets, the thermal axion mass 
constraints are degraded only slightly. The addition of the measurements of as and Q m from the 
2013 Planck cluster catalogue on galaxy number counts relaxes the bounds on the thermal axion 
mass, mildly favouring a ~ 1 eV axion mass, regardless of the model adopted for the primordial 
power spectrum. However, in general, such a preference disappears if the sum of the three active 
neutrino masses is also considered as a free parameter in our numerical analyses, due to the strong 
correlation between the masses of these two hot thermal relics. 
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I. INTRODUCTION 


A possible candidate for an extra hot thermal relic 
component is the axion particle produced thermally in 
the early universe. Axions therefore can contribute to 
the hot dark matter component together with the stan¬ 
dard relic neutrino background. Axions may be produced 
in the early universe via thermal or non thermal pro¬ 
cesses, and arise as the solution to solve the strong CP 
problem w- Axions are the Pseudo-Nambu-Goldstone 
bosons of a new global U( 1) pq (Peccei-Quinn) symmetry 
that is spontaneously broken at an energy scale f a . The 
axion mass is given by 


/ttTOtt y/R „ „ „ 10 7 GeV 
m “ = — TTR=°- 6eV -^ 


(i) 


where f a is the axion coupling constant, R = 0.553 ± 
0.043 is the up-to-down quark masses ratio and f n = 93 
MeV is the pion decay constant. Non-thermal axions, as 
those produced by the misalignment mechanism, while 
being a negligible hot dark matter candidate, may con¬ 
stitute a fraction or the total cold dark matter compo¬ 
nent of the universe. We do not explore such a possi¬ 
bility here. Thermal axions will affect the cosmological 
observables in a very similar way to that induced by the 
presence of neutrino masses and/or extra sterile neutrino 
species. Massive thermal axions as hot relics affect large 
scale structure, since they will only cluster at scales larger 
than their free-streaming scale when they become non- 
relativistic, suppressing therefore structure formation at 
small scales (large wavenumbers k). Concerning Cosmic 
Microwave Background (CMB) physics, an axion mass 
will also lead to a signature in the CMB photon temper¬ 
ature anisotropies via the early integrated Sachs-Wolfe 
effect. In addition, extra light species as thermal axions 


will contribute to the dark radiation content of the uni¬ 
verse, or, in other words, will lead to an increase of the 
effective number of relativistic degrees of freedom N e g, 
defined via 


Prad — 




P"1 ? 


( 2 ) 


where p 7 refers to the present photon energy density. 
In the standard cosmological model in which a thermal 
axion content is absent, the three active neutrino con¬ 
tribution leads to the canonical value of N e s = 3.046 
dj. The extra contribution to N e g arising from thermal 
axions can modify both the CMB anisotropies (via Silk 
damping) and the light element primordial abundances 
predicted by Big Bang Nucleosynthesis. The former cos¬ 
mological signatures of thermal axions have been exten¬ 
sively exploited in the literature to derive bounds on the 
thermal axion mass, see Refs. mm- 

However, all the cosmological axion mass limits to date 
have assumed the usual simple power-law description for 
the primordial perturbations. The aim of this paper is 
to constrain the mass of the thermal axion using a non- 
parametric description of the Primordial Power Spectrum 
(PPS hereinafter) of the scalar perturbations, as intro¬ 
duced in Ref. m- While in the simplest models of in¬ 
flation jl2fi22] the PPS has a scale-free power-law form, 
the PPS could be more complicated, presenting various 
features or a scale dependence. Several methods have 
been proposed in the literature to reconstruct the shape 
of the PPS (see the recent work of Ref. 23]). It has 
been shown (23[ E2>| that there are small hints for devia¬ 
tions from the power-law form, even when using different 
methods and different data sets. 

The energy scales at which the PPS was produced dur¬ 
ing inflation can not be directly tested. We can only infer 






2 


the PPS by measuring the current matter power spec¬ 
trum in the galaxy distribution and the power spectrum 
of the CMB fluctuations. The latter one, measured with 
exquisite precision by the Planck experiment |26H28| . is 
the convolution of the PPS with the transfer function. 
Therefore, in order to reconstruct the PPS, the assump¬ 
tion of an underlying cosmological model is a mandatory 
first step in order to compute the transfer function. 

Here we rather exploit a non-standard PPS approach, 
which can allow for a good fit to experimental data even 
in models that deviates from the standard cosmological 
picture. In particular, we consider a thermal axion sce¬ 
nario, allowing the PPS to assume a more general shape 
than the usual power-law description. This will allow us 
to test the robustness of the cosmological thermal axion 
mass bounds (see Ref. [ID] for a recent standard thermal 
axion analysis), as first performed in Ref. |29j for the 
neutrino mass case. 

The structure of the paper is as follows. Section |TT] de¬ 
scribes the PPS modeling used in this study, as well as the 
description of the thermal axion model explored here and 
the cosmological data sets exploited to constrain such a 
model. In Sec. [nT]we present and discuss the results aris¬ 
ing from our bayesian analysis, performed through the 
Monte Carlo Markov Chains (MCMC) package CosmoMC 
ED], while the calculation of the theoretical observables is 
done through the Boltzman equations solver CAMB (Code 
for Anisotropies in the Microwave Background) [31]. We 
draw our conclusions in Sec. nvi 


tion to interpolate the PPS values in a series of nodes 
at fixed position. The interpolating function we used 
is named “piecewise cubic Hermite interpolating polyno¬ 
mial” (PCHIP) [33] and it is a modified spline function, 
defined to preserve the original monotonicity of the point 
series that is interpolated. We use a modified version of 
the original PCHIP algorithm [35] , detailed in Appendix A 
of Ref. [XT . 

To describe the scalar PPS with the PCHIP model, we 
only need to give the values of the PPS in a discrete 
number of nodes and to interpolate among them. We 
use 12 nodes which span a wide range of k- values: 

ki= 5 x 10 -6 Mpc -1 , 

k '2 = 10 - ' 1 Mpc -1 , 

kj = fc 2 (fcii/fc 2 ) (j-2)/9 for j G [3,10], 
fcn = 0.35 Mpc -1 , 

ki 2 = 10 Mpc -1 . (3) 

We choose equally spaced nodes in the logarithmic scale 
in the range (/c 2 ,fcn), that is well constrained from the 
data [2D], while the first and the last nodes are useful to 
allow for a non-constant behaviour of the PPS outside 
the well-constrained range. 

The PCHIP PPS is described by 

P s (k) =P 0 x PCHIP(fc; P 8|1 ,..., P Sl i 2 ), (4) 

where P s j is the value of the PPS at the node kj divided 
by P 0 = 2.36 x 10 -9 [SB]. 


II. METHOD 


B. Cosmological and Axion Model 


In this section we focus on the tools used in the numer¬ 
ical analyses performed here. Subsection II A| describes 
the alternative model for the PPS of scalar perturbations 
used for the analyses here (see also Ref. [13]), while in 
Subsection |II B| we introduce the cosmological model and 
the thermal axion treatment followed in this study. Fi¬ 
nally, we shall present in Subsection |II C| the cosmological 
data sets used in the MCMC analyses. 


A. Primordial Power Spectrum Model 

The primordial fluctuations in scalar and tensor modes 
are generated during the inflationary phase in the early 
universe. The simplest models of inflation predict a 
power-law form for the PPS of scalar and tensor per¬ 
turbations (see e.g. [T2H32| and references therein), but 
in principle inflation can be generated by more compli¬ 
cated mechanisms, thus giving a different shape for the 
PPS (see Refs. [32] S3] and references therein). In order 
to study how the cosmological constraints on the param¬ 
eters change in more general inflationary scenarios, we 
assume a non-parametric form for the PPS. 

Among the large number of possibilities, we decided 
to describe the PPS of scalar perturbations using a func- 


The baseline scenario we consider here is the ACDM 
model, extended with hot thermal relics (the axions), to¬ 
gether with the PPS approach outlined in the previous 
section. For the numerical analyses we use the following 
set of parameters, for which we assume flat priors in the 
intervals listed in Tab. [3 

{ujb,w c , @ s ,t, . ,P a , 12 } , (5) 

where = flbh 2 and oj c = Vt c h 2 are, respectively, the 
physical baryon and cold dark matter energy densities, 
©s is the ratio between the sound horizon and the an¬ 
gular diameter distance at decoupling, r is the reioniza¬ 
tion optical depth, m a and T) m v are the axion and the 
sum of three active neutrino masses (both in eV) and 
P s i,..., P a i 2 are the parameters of the PCHIP PPS. We 
shall also consider a scenario in which massive neutrinos 
are also present, to explore the expected degeneracy be¬ 
tween the sum of the neutrino masses and the thermal 
axion mass m- 

In order to compare the results obtained with the 
PCHIP PPS to the results obtained with the usual power- 
law PPS model, we describe the latter case with the fol¬ 
lowing set of parameters: 

cu c , ©s, r, 7?z a , ti s , log[10 10 A s ]} , (6) 
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Parameter 

Prior 

fib h 2 

[0.005,0.1] 

^cdm^ 

[0.001,0.99] 

e s 

[0.5,10] 

T 

[0.01,0.8] 

m a (eV) 

[0.1,3] 

( eV ) 

[0.06,3] 

Pa, 1, ■ • • , Ps, 12 

[0.01,10] 

n s 

[0.9,1.1] 

log[10 10 A s ] 

[2.7,4] 


TABLE I: Priors for the parameters used in the MCMC anal¬ 
yses. 


where n s is the scalar spectral index, A s the amplitude 
of the primordial spectrum and the other parameters are 
the same ones described above. The case of several hot 
thermal relics for the standard scenario will not be car¬ 
ried out here, as it has been done in the past by several 
authors (see e.g. [IQ]). The flat priors we use are listed 
in Tab. [U 

Concerning the contribution of the axion mass-energy 
density to the universe’s expansion rate, we briefly sum¬ 
marize our treatment in the following. Axions decoupled 
in the early universe at a temperature Tjj given by the 
usual freeze out condition for a thermal relic: 


T(T C ) = H(T d ) , 


(7) 


where the thermally averaged interaction rate T refers to 
the 7 r T 7 r —> 7 r + a process: 


T = 


1 


10247T 5 flfl 


r 2 t 

( - / an 1 a 


2 t 2 a7r 


( 8 ) 


with C a -K = 3 ( 1 i + J r) representing the axion-pion coupling 
constant and the integral I a reads as follows 


„2„,2 


I a = n a T s / dxidx 2 - 




y W 2 

I" 1 


s 2 T 4 


(9) 


in which n a = (£ 3 / 7 r 2 )T 3 refers to the number density 
for axions in thermal equilibrium. The function f(y) = 
l/(e y — 1 ) is the pion thermal distribution and there are 
three different kinematical variables (xi = \pi\/T, yi = 
Ei/T (i = 1,2) and s = 2(ml + T 2 (y 1 y 2 -XiX 2 oj))). The 
freeze out equation above, Eq. |7]), can be numerically 
solved, obtaining the axion decoupling temperature Tjj 
as a function of the axion mass m a . Figure [l] shows, in 
the left upper panel, the axion decoupling temperature 
as a function of the axion mass, in eV units. Notice that, 
the higher the axion mass, the lower the temperature of 
decoupling is. From the axion decoupling temperature 
it is possible to infer the present axion number density, 
related to the current photon density ra 7 by 


_ g*s(To ) n 7 
Ua ~ 9 *s(T d ) X 2 


( 10 ) 


Axion parameter 



m a (eV) 

0.1 

3 

T d (MeV) 

245.6 

43.2 

Q a h 2 

0.0003 

0.016 

AiVeff 

0.18 

0.43 

fcf s (/i/Mpc) 

0.06 

1.46 


TABLE II: Values for the axion parameters, To, AN e g, fl a h 2 
and fcfs for the lower and upper prior choice of m a explored 
here. 


where g+s represents the number of entropic degrees of 
freedom, with g*s(To) = 3.91. The contribution of the 
relic axion to the total mass-energy density of the uni¬ 
verse will be given by the product of the axion mass 
times the axion number density. The quantity tt a h 2 at 
the present epoch is depicted in the bottom left panel of 
Fig- [l] Notice that, currently, a 1 eV axion will give rise 
to f i a h 2 — 0.005, while a neutrino of the same mass will 
contribute to the total mass-energy density of the uni¬ 
verse with f l u h 2 ~ 0.01. Notice however that f 1 a h 2 rep¬ 
resents the contribution from relic, thermal axion states. 
Non-thermal processes, as the misalignment production, 
could also produce a non-thermal axion population which 
we do not consider here, see the work of m for the most 
recent cosmological constraints on such scenario. As pre¬ 
viously stated, the presence of a thermal axion will also 
imply an extra radiation component at the BBN period: 


AlVeJ = 


4 

7 


/ 3 rioA 4/3 
\2nJ 


( 11 ) 


where n a is given by Eq. ( |10[ ) and n v refers to the present 
neutrino plus antineutrino number density per flavour. 
The top right panel of Fig. [T] shows the axion contri¬ 
bution to the radiation component of the universe as a 
function of the axion mass. Notice that the extra dark 
radiation arising from a 1 eV axion is still compatible 
(at 95% CL) with the most recent measurements of N e fi 
from the Planck mission [26]. The last crucial cosmo¬ 
logical axion quantity is the axion free streaming scale, 
i.e. the wavenumber k{ s below which axion density per¬ 
turbations will contribute to clustering once the axion 
is a non-relativistic particle. This scale is illustrated in 
Fig.0 in the bottom right panel, together with that cor¬ 
responding to a neutrino of the same mass. Notice that 
they cover the same scales for our choice of priors for m a 
and an< l therefore one can expect a large corre¬ 

lation between these two quantities in measurements of 
galaxy clustering. We will explore this degeneracy in the 
following sections. We summarize the axion parameters 
in Tab. [TTJ where we specify the values of the decoupling 
temperature, A7V e fj, fl a h 2 and kf s for the range of axion 
masses considered here, (0.1,3) eV. 
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FIG. 1: The left upper panel shows the temperature of decoupling as a function of the axion mass (solid curve), as well as the 
Big Bang Nucleosynthesis temperature, Tbbn — 1 MeV (dashed curve). The right upper panel shows the axion contribution to 
the extra dark radiation content of the universe, while the bottom right plot depicts the free-streaming scale of an axion (solid 
curve) or a neutrino (dashed curve) versus the axion/neutrino mass, in eV. The left bottom panel shows the current axion 
mass-energy density as a function of the axion mass. 


C. Cosmological measurements 

Our baseline data set consists of CMB measurements. 
These include the temperature data from the Planck 
satellite, see Refs. [MISS], together with the WMAP 9- 
year polarization measurements, following [39] . We also 
consider high multipole data from the South Pole Tele¬ 
scope (SPT) [3(]j] as well as from the Atacama Cosmology 
Telescope (ACT) [41] . The combination of all the above 
CMB data is referred to as the CMB data set. 

Galaxy clusters represent an independent tool to probe 
the cosmological parameters. Cluster surveys usually re¬ 
port their measurements by means of the so-called clus¬ 
ter normalization condition, where 7 ~ 0.4 [424 

[33] , We shall use here the cluster normalization condi¬ 
tion as measured by the Planck Sunyaev-Zeldovicli (PSZ) 
2013 catalogue |4£], referring to it as the PSZ data set. 
The PSZ measurements of the cluster mass function pro¬ 
vide the constraint cr$(£l m / 0.27) 0 3 = 0.764 ± 0.025. As 
there exists a strong degeneracy between the value of 
the ag parameter and the cluster mass bias, it is possi¬ 
ble to fix the value of the bias parameter accordingly to 
the results arising from numerical simulations. In this 
last case, the error on the cluster normalization con¬ 


dition from the PSZ catalogue is considerably reduced: 
O8(^m/0.27) 0 ’ 3 = 0.78 ± 0.01. In our analyses, we shall 
consider the two PSZ measurements of the cluster nor¬ 
malization condition, to illustrate the impact of the clus¬ 
ter mass bias in the thermal axion mass bounds, as re¬ 
cently explored in Ref. j46j for the neutrino mass case. 
Figure [2] illustrates the prediction for the cluster normal¬ 
isation condition, as(Clm/ 0.27) 03 , as a function of the 
thermal axion mass. We also show the current PSZ mea¬ 
surements with their associated 95% CL uncertainties, in¬ 
cluding those in which the cluster mass bias parameter is 
fixed. Notice that the normalisation condition decreases 
as the axion mass increases, due to the decrease induced 
in the a$ parameter in the presence of axion masses: the 
larger the axion mass is, the larger the reduction in the 
matter power spectra will be. 

Tomographic weak lensing surveys are sensitive to the 
overall amplitude of the matter power spectrum by mea¬ 
suring the correlations in the observed shape of distant 
galaxies induced by the intervening large scale struc¬ 
ture. The matter power spectrum amplitude depends on 
both the as clustering parameter and the matter den¬ 
sity fl m . Consequently, tomographic lensing surveys, 
via measurements of the galaxy power shear spectra, 
provide additional and independent constraints in the 
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m a [eV] 


FIG. 2: Cluster normalisation condition, crg(n m /0.27)°' 3 , as a function of the thermal axion mass. We also show the current 
PSZ measurements with their associated 95% CL uncertainties. 


(as, H m ) plane. The Canada-France-Hawaii Telescope 
Lensing Survey, CFHTLenS, with six tomographic red- 
shift bins (from 2 = 0.28 to z = 1.12), provides a 
constraint on the relationship between ag and il m of 
cr 8 (fl m /0.27)°- 46 = 0.774 ± 0.040 g7]. We shall refer to 
this data set as CFHT. 

We also address here the impact of a gaussian prior on 
the Hubble constant Hq = 70.6 ±3.3 km/s/Mpc from an 
independent reanalysis of Cepheid data [48] . referring to 
this prior as the HST data set. 

We have also included measurements of the large scale 
structure of the universe in their geometrical form, i.e., 
in the form of Baryon Acoustic Oscillations (BAO). Al¬ 
though previous studies in the literature have shown that, 
for constraining hot thermal relics, the shape information 
contained in the galaxy power spectrum is more power¬ 
ful when dealing simultaneously with extra relativistic 
species and hot thermal relic masses S3 Eg, we exploit 
here the BAO signature, as the contribution from the 
thermal axions to the relativistic number of species is 
not very large (see Tab. |H|, and current measurements 
from galaxy surveys are mostly reported in the geomet¬ 
rical (BAO) form. 

The BAO wiggles, imprinted in the power spectrum of 
the galaxy distribution, result from the competition in 
the coupled photon-baryon fluid between radiation pres¬ 
sure and gravity. The BAO measurements that have been 
considered in our numerical analyses include the results 
from the WiggleZ [5T|, the 6dF [52] and the SDSS II sur¬ 
veys [55] [54], at redshifts of z = 0.44,0.6,0.73, z = 0.106 
and 2 = 0.35, respectively. We also include in our analy¬ 
ses as well the Data Release 11 (DR11) BAO signal of the 
BOSS experiment [55], which provides the most precise 
distant constraints [55 measuring both the Hubble pa¬ 
rameter and the angular diameter distance at an effective 


redshift of 0.57. Figure [3] illustrates the spherically aver¬ 
aged BAO distance, Dy(z) oc D\(z)/H(z), as a function 
of the axion mass, at a redshift of z = 0.57, as well as the 
measurement from the BOSS experiment with 95% CL 
error bars [56] . Notice that, from background measure¬ 
ments only, there exists a strong degeneracy between the 
cold dark matter mass-energy density and the axion one. 
The solid black line in Fig. [3] shows the spherically av¬ 
eraged BAO distance if all the cosmological parameters 
are fixed, including u> c . The spherically averaged BAO 
distance deviates strongly from the ACDM prediction. 
However, if uj c is varied while m a is changed (in order 
to keep the total matter mass-energy density constant, 
see the dotted blue line in Fig. [3]), the spherically aver¬ 
aged BAO distance approaches to its expected value in a 
ACDM cosmology. 

D. Compatibility of data 

It has been pointed out (see Sec. 5.5 of Ref. m and 
also Refs. HUSH) that the value of erg reported by clus¬ 
ter measurements and the value estimated from Planck 
CMB measurements show a tension at the ~ 2cr level. 
These discrepancies may arise due to the lack of a full 
understanding of the cluster mass calibrations. Although 
some studies in the literature, including the present one, 
show that in extended cosmological models with non-zero 
neutrino masses the discrepancies previously mentioned 
could be alleviated, the results from Ref. m show, us¬ 
ing also Bayesian evidence, that a canonical ACDM sce¬ 
nario with no massive neutrinos is preferred over its neu¬ 
trino extensions by several combinations of cosmological 
datasets. Therefore, the results presented here and ob¬ 
tained when considering cluster data depend strongly on 
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FIG. 3: The solid black line depicts the spherically averaged BAO distance Dy{z), as a function of the axion mass, at a redshift 
of z = 0.57, after keeping fixed all the remaining cosmological parameters, the cold dark matter included. The dashed blue 
line depicts the equivalent but keeping fixed the total matter mass-energy density (and consequently changing the cold dark 
matter w c ). The bands show the measurement from the BOSS experiment (DR11) with its associated 95% CL error. 


the reliability of low-redshift cluster data. If future data 
confirm current low-redshift cluster measurements, one 
could further test some of the possible beyond the ACDM 
models using particle physics experiments. For instance, 
the existence of a full thermal sterile neutrino could be 
tested with neutrino oscillation experiments, and the ac¬ 
tive neutrino mass could also be tested by tritium ex¬ 
periments or, if the neutrino is a Majorana particle, by 
neutrinoless double beta decay searches. 


III. RESULTS 


Table III depicts our results in the first scenario ex¬ 


plored here, in which the axion mass is a free parame¬ 
ter and the PPS is described by the approach specified 
in Sec. |II A| Concerning CMB measurements only, the 
bounds on the thermal axion masses are largely relaxed 
in the case in which the PPS is not described by a simple 
power-law, as can be noticed after comparing the results 
depicted in Tab. |III| with those shown in Tab. |IV| This 
can be understood in terms of Fig. [4j which illustrates the 
degeneracies in the temperature anisotropies between the 
thermal axion mass and the PCHIP PPS. Figure [4] shows 
the temperature anisotropies for a ACDM model and a 
power-law PPS (solid red line), for a 2 eV thermal axion 
mass and a power-law PPS (dashed blue line) and for a 
ACDM model but the PPS described by the PCHIP model 
explored here (dotted black line), with values for the P s j 
chosen to match the non-zero thermal axion mass curve, 


accordingly to their current allowed regions (see Tab. Ill I. 


More concretely, we have used the following values for the 
PPS parameters: P s q = 1.15, P Sy 2 = 1.073, P s ,3 = 1.058, 
Ps ,4 = 1-03, P s , 5 = 0.99, P afi ’= 0.97, pj = 0.966, 
P Sj8 = 0.932, P ,, 9 = 0.91, P a ,i 0 = 0.86, P s> n = 0.84 and 
P, 12 = 0.77. Notice that the case of a 2 eV thermal ax¬ 
ion can be easily mimicked by a simple ACDM model if 


the assumptions concerning the PPS shape are relaxed. 
We also add in this figure the measurements of the pho¬ 
ton temperature anisotropies from the Planck 2013 data 
release m 

The addition to the CMB data of the HST prior on 
the Hubble constant provides a 95% CL upper limit on 
the thermal axion mass of 1.31 eV 1 , while the further 
addition of the BAO measurements brings this constraint 
down to 0.91 eV, as these last data sets are directly sen¬ 
sitive to the free-streaming nature of the thermal axion. 
Notice that these two 95% CL upper bounds are very 
similar to the ones obtained when considering the stan¬ 
dard power-law power spectrum, which are 1.56 eV and 
0.83 eV for the CMB+HST and CMB+HST+BAO data 
combinations, respectively. 

Interestingly, when adding the CFHT bounds on the 
<r 8 -fl m relationship, the bounds on the thermal axion 
mass become weaker. The reason for that is due to the 
lower <r 8 values preferred by weak lensing measurements, 
values that can be achieved by allowing for higher ax¬ 
ion masses. The larger the axion mass, the larger is the 
reduction of the matter power spectrum at small (i.e. 
cluster) scales, leading consequently to a smaller value of 
the clustering parameter cr 8 . 

If we instead consider now the PSZ data set with fixed 
cluster mass bias, together with the CMB, BAO and HST 
measurements, a non-zero value of the thermal axion 
mass of ~ 1 eV (~ 0.80 eV) is favoured at ~ 4er (~ 3a) 
level, when considering the PCHIP (standard power-law) 


1 There exists a very large degeneracy between H 0 and the neu¬ 
trino masses when restricting the numerical analyses to CMB 
measurements. The addition of the HST prior on the Hubble 
constant helps enormously in breaking this degeneracy, see W- 
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TABLE III: 95% CL constraints on the physical cold dark matter density Q c h 2 , the axion mass m a (in eV), the clustering 
parameter og, the relative matter energy density f2 m and the P s j parameters for the PPS nodes from the different combinations 
of data sets explored here in the ACDM+m a model, considering the PCHIP PPS modeling. 


PPS approach 2 . However, these results must be re¬ 
garded as an illustration of what could be achieved with 
future cluster mass calibrations, as the Planck collabo¬ 
ration has recently shown in their analyses of the 2015 
Planck cluster catalogue [2B]. When more realistic ap¬ 
proaches for the cluster mass bias are used, the errors on 
the so-called cluster normalization condition are larger, 
and, consequently, the preference for a non-zero axion 
mass of 1 eV is only mild in the PCHIP PPS case, while 
in the case of a standard power-law PPS such an evidence 
completely disappears. 

Figure [5] (left panel) shows the 68% and 95% CL al¬ 
lowed regions in the (m a , H c /i 2 ) plane for some of the 
possible data combinations explored in this study, and 
assuming the PCHIP PPS modeling. Notice that, when 
adding BAO measurements, lower values of the physical 


2 


A similar effect when considering PSZ data for constraining ei¬ 
ther thermal axion or neutrino masses has also been found in 


Refs. Hol I58H6T1 . 


cold dark matter density are preferred. This is due to 
the fact that large scale structure allows for lower axion 
masses than CMB data alone. The lower is the thermal 
axion mass, the lower is the amount of hot dark matter 
and consequently the lower should be the cold dark mat¬ 
ter component. This effect is clear from the results shown 


in Tab. Ill and Tab. |IV| where the values of the physical 
cold dark matter density f l c h 2 and of the relative cur¬ 
rent matter density arising from our numerical fits 
are shown, for the different data combinations considered 
here. 

The right panel of Fig. [5] shows the 68% and 95% CL 
allowed regions in the (m Q , erg) plane in the PCHIP PPS 
scenario. The lower values of the erg clustering parameter 
preferred by PSZ data (see the results shown in Tab. Ill 
and Tab. IVI are translated into a preference for non-zero 


thermal axion masses. Larger values of m a will enhance 
the matter power spectrum suppression at scales below 
the axion free-streaming scale, leading to smaller values 
of the (jg clustering parameter, as preferred by PSZ mea¬ 
surements. The evidence for non-zero axion masses is 
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TABLE IV: 95% CL constraints on f 1 c h 2 , the axion mass m a (in eV), og, fi m , log(10 1(, A s ) and n s from the different combinations 
of data sets explored here in the ACDM+m a model, assuming the standard power-law PPS. 
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FIG. 4: Temperature anisotropies for the pure ACDM model and a power-law PPS (solid red line), for a 2 eV thermal axion 
mass and a power-law PPS (dashed blue line) and for the standard ACDM model but the PPS described by the PCHIP model 
(dotted black line). The data points and the error bars in the left panel show the measurements of the photon temperature 
anisotropies arising from the Planck 2013 data release 123 - 


more significant when fixing the cluster mass bias in the 
PSZ data analyses. 

Figure [6] shows the equivalent to Fig. [5] but for a stan¬ 
dard power-law PPS. Notice that, except for the case in 
which CMB measurements are considered alone, the ther¬ 
mal axion mass constraints do not change significantly, if 
they are compared to the PCHIP PPS modeling. This fact 
clearly states the robustness of the cosmological bounds 
on thermal axion masses and it is applicable to the re¬ 
maining cosmological parameters, see Tabs. |III| and |IV| 
Note that, for the standard case of a power-law PPS, 
the preference for non-zero axion masses appears only 
when considering the (unrealistic) PSZ analysis with a 
fixed cluster mass bias. When more realistic PSZ mea¬ 
surements of the cluster normalization condition are ex¬ 
ploited, there is no preference for a non-zero thermal ax¬ 


ion mass. 

The last scenario we explore here is a 
ACDM+m a +^ universe, in which we consider 
two coexisting hot dark matter species: thermal axions 
and three active (massive) neutrinos. Table [v| illustrates 
the equivalent of Tab. |III| but including the active 
neutrino masses in the MCMC parameters. We do 
not perform here the analysis for the hot mixed dark 
matter model with the standard power-law matter 
power spectrum, as it was already presented previously 
in Ref. [TO] . If we compare to the standard power-law 
case, we find that the bounds on the axion and neutrino 
masses presented here are very similar. Furthermore, no 
evidence for neutrino masses nor for a non-zero axion 
mass appears in this mixed hot dark matter scenario 
(except for the axion case and only if considering PSZ 
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FIG. 5: The left panel depicts the 68% and 95% CL allowed regions in the ( m a , £l c h 2 ) plane for different possible data 
combinations, when a PCHIP PPS is assumed. The right panel shows the equivalent but in the ( m a , erg) plane. 
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FIG. 6: The left panel depicts the 68% and 95% CL allowed regions in the (m Q , fl c h 2 ) plane for different possible data 
combinations, when a power-law PPS is assumed. The right panel shows the equivalent but in the ( m a , as) plane. 


clusters with the bias fixed). The reason for that is 
due to the strong degeneracy between m a and 
see Fig. [7J where one can notice that that these two 
parameters are negatively correlated: an increase in the 
axion mass will increase the amount of the hot dark 
matter component. In order to compensate the changes 
in both the CMB temperature anisotropies (via the 
early ISW effect) and in the power spectrum (via the 


suppression at small scales of galaxy clustering), the 
contribution to the hot dark matter from the neutrinos 
should be reduced. We have shown in Fig. [7] three 
possible data combinations. Notice that for the case in 
which PSZ cluster measurements (with the bias fixed) 
are included the strong degeneracy between m a and 
m„ is partially broken, due to the smaller value of ag 
preferred by the former data set. However, these results 
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FIG. 7: 68% and 95% CL allowed regions in the (X) m "i m “) 
plane, both in eV, for three different possible data combina¬ 
tions, when a PCHIP PPS is assumed. 


strongly rely on the numerical results concerning the 
cluster mass bias and therefore the evidence for m a ^ 0 
should be regarded as what could be obtained in if these 
measurements are further supported by independent 
data from future cluster surveys. 


Besides the results concerning the thermal axion mass 
and the standard ACDM parameters, we also obtain con¬ 
straints on the form of the PPS when modeled accord¬ 
ingly to the PCHIP scenario. The 95% CL limits for the 
P s :j parameters are shown in Tab. Ill while an exam¬ 
ple of the reconstructed PPS is given in Fig. [HJ where 
we show the 68%, 95% and 99% CL allowed regions aris¬ 
ing from a fit to CMB data of the PCHIP PPS scale de¬ 
pendence, in the context of a ACDM+m a model. We 
do not show the corresponding figures obtained from all 
the other data combinations since they are equivalent to 
Fig-0 as one can infer from the very small differences in 
the 95% CL allowed ranges for the P s j parameters aris¬ 
ing from different data sets, see Tab. Ill Note that both 
P Sj i and -P s , 12 are poorly constrained at this confidence 
level: the reason for that is the absence of measurements 
at their corresponding wavenumbers. All the remaining 
P s j. with j = 2,..., 11 are well-constrained. In par¬ 
ticular, in the range between fc 5 and k w (see Eq. (|3|), 
the P S J are determined with few percent accuracy. In¬ 
deed, in the range covered between these nodes, the 
PPS does not present features and can be perfectly de¬ 
scribed by a power-law parametrization. Among the in¬ 
teresting features outside the former range, we can no¬ 
tice in Fig. [8] a significant dip at wavenumbers around 


k = 0.002 Mpc -1 , that comes from the dip at i = 20 — 30 
in the CMB temperature power spectrum and a small 
bump around k = 0.0035 Mpc -1 , corresponding to the 
increase at £ ~ 40. These features have been obtained in 
previous works OEUS using different methods and 
data sets. In addition, we obtain an increase of power at 
k ~ 0.2 Mpc -1 , necessary to compensate the effects of the 
thermal axion mass in both the temperature anisotropies 
and the large scale structure of the universe. 


IV. CONCLUSIONS 

Axions provide the most elegant scenario to solve the 
strong CP problem, and may be produced in the early 
universe via both thermal and non-thermal processes. 
While non thermal axions are highly promising cold dark 
matter candidates, their thermal companions will con¬ 
tribute to the hot dark matter component of the universe, 
together with the (light) three active neutrinos of the 
standard model of elementary particles. Therefore, the 
cosmological consequences of light massive thermal ax¬ 
ions are very much alike those associated with neutrinos, 
as axions also have a free-streaming nature, suppressing 
structure formation at small scales. Furthermore, these 
light thermal axions will also contribute to the dark radi¬ 
ation background, leading to deviations of the relativis¬ 
tic degrees of freedom N e g from its canonically expected 
value of N e g = 3.046. Based on these signatures, several 
studies have been carried out in the literature deriving 
bounds on the thermal axion mass mm- 

Nevertheless, these previous constraints assumed that 
the underlying primordial perturbation power spectrum 
follows the usual power-law description governed, in its 
most economical form, by an amplitude and a scalar spec¬ 
tral index. Here we have relaxed such an assumption, in 
order to test the robustness of the cosmological axion 
mass bounds. Using an alternative, non-parametric de¬ 
scription of the primordial power spectrum of the scalar 
perturbations, named PCHIP and introduced in Ref. mi. 
we have shown that, in practice, when combining CMB 
measurements with low redslrift cosmological probes, the 
axion mass constraints are only mildly sensitive to the 
primordial power spectrum choice and therefore are not 
strongly dependent on the particular details of the un¬ 
derlying inflationary model. These results agree with the 
findings of Ref. [25] for the neutrino mass case. The tight¬ 
est bound we find in the PCHIP primordial power spec¬ 
trum approach is obtained when considering BAO mea¬ 
surements together with CMB data, with m a < 0.89 eV 
at 95% CL. In the standard power-law primordial power 
spectrum modeling, the tightest bound is m a < 0.83 eV 
at 95% CL, obtained when combining BAO, CMB and 
HST measurements. Notice that these bounds are very 
similar, confirming the robustness of the cosmological ax¬ 
ion mass measurements versus the primordial power spec¬ 
trum modeling. 

Interestingly, both weak lensing measurements and 
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TABLE V: 95% CL constraints on the physical cold dark matter density Q c h 2 , the axion mass m a , the sum of the active neutrino 
masses (both in eV), the clustering parameter as, the relative matter energy density fi m and the P s j parameters for the 

PPS nodes from the different combinations of data sets explored here in the ACDM+m a +^ rth, model, considering the PCHIP 
PPS modeling. 


cluster number counts weaken the thermal axion mass 
bounds. The reason for that is due to the lower as values 
preferred by these measurements, which could be gener¬ 
ated by a larger axion mass. More concretely, Planck 
cluster measurements provide a measurement of the so- 
called cluster normalization condition, which establishes 
a relationship between the clustering parameter cr 8 and 
the current matter mass-energy density f2 m . However, 
the errors on this relationship depend crucially on the 
knowledge of the cluster mass bias. A conservative ap¬ 
proach for the cluster mass calibration results in a mild 
(zero) evidence for a non-zero axion mass of 1 eV in 
the PCHIP (power-law) PPS case. We also illustrate a 
case in which the cluster mass bias is fixed, to forecast 
the expected results from future cosmological measure¬ 
ments. In this case, a non-zero value of the thermal ax¬ 
ion mass of ~ 1 eV (~ 0.80 eV) is favoured at ~ 4cr 
(~ 3tr) level, when considering the PCHIP (power-law) 
PPS approach. When considering additional hot relics 


in our analyses, as the sum of the three active neutrino 
masses, the evidence for a ~ 1 eV thermal axion mass 
disappears almost completely. Furthermore, these values 
of axion masses correspond to an axion coupling constant 
f a = 6 x 10 6 GeV, which seems to be in tension with the 
limits extracted from the neutrino signal duration from 
SN 1987A [B51|53] (albeit these limits depend strongly on 
the precise axion emission rate and still remain rough es¬ 
timates). Precise cluster mass calibration measurements 
are therefore mandatory to assess whether there exists a 
cosmological indication for non-zero axion masses, as the 
cluster mass bias is highly correlated with the clustering 
parameter erg, which, in turn, is highly affected by the 
free-streaming nature of a hot dark matter component, 
as thermal axions. 
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FIG. 8: 68%, 95% and 99% CL allowed regions for the PCHIP PPS scale dependence in the ACDM+m a model, using CMB data 
only. The bands are obtained with a marginalization of the posterior distribution for each different value of the wavenumber k 
in a fine grid. The black line represents the peak of the posterior distribution at each value of k. 
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